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Networks of globally coupled, noise activated, bistable elements with connection time delays are 
considered. The dynamics of these systems is studied numerically using a Langevin description 
and analytically using (1) a Gaussian approximation as well as (2) a dichotomous model. The 
system demonstrates ordering phase transitions and multi-stability. That is, for a strong enough 
feedback it exhibits nontrivial stationary states and oscillatory states whose frequencies depend only 
on the mean of the time delay distribution function. Other observed dynamical phenomena include 
coherence resonance and, in the case of non-uniform coupling strengths, amplitude death and chaos. 
Furthermore, an increase of the stability of the trivial equilibrium with increasing non-uniformity 
of the time delays is observed. 
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I. INTRODUCTION 

Due to its relevance for a variety of scientific disciplines 
such as physics, chemistry, biology, economics and social 
sciences, the study of collective phenomena in extended 
stochastic systems with long range interaction has been 
of great interest in recent years and various techniques 
based on Langevin, Fokker-Planck and master equations 
have been conceived to explore their dynamics. 

An effective and simple model for the study of noise 
induced collective phenomena is the globally coupled net- 
work of stochastically driven bistable elements. Indeed, 
the cooperative dynamics of these systems has been the 
subject of many studies and its relevance for critical phe- 
nomena 0, spin systems H, neural networks [3|.[J[5|. 
genetic regulatory networks [6j and decision making pro- 
cesses in social systems Q has been pointed out. 

For the sake of simplicity it has traditionally been as- 
sumed that the interactions in these networks are instan- 
taneous. However, in recent years it has been realized 
that time delays due to finite transmission and process- 
ing speeds are 1) significant compared to the dynamical 
time scales of the system and 2) often change fundamen- 
tally its dynamical properties H, 0, 0, . 

Thus, in this paper the generic model of globally cou- 
pled bistable elements is extended by time delayed cou- 
plings and its collective dynamics is studied numerically 
and analytically. 

The properties of globally coupled dynamical units, rel- 
evant for system such as arrays of lasers 0] an d Joseph- 
son junctions 1131, h ave been explored in many studies 
[see alsollQElMiniEl- Desai and Zwanzig [3, for 
instance, studied the synchronization of noise activated 
bistable oscillators with instantaneous coupling and de- 
rived from the Fokker-Planck equation for the joint prob- 
ability distribution of the oscillators an exact mean field 
model (DZ-model) in the thermodynamic limit N — > oo, 
where N is the number of network elements. Beyond a 
critical coupling strength this system displays a second 
order phase transition to an ordered nontrivial station- 
ary state. The effect of uniform interaction delays in a 



globally coupled network of phase oscillators has been 
explored by Yeung and Strogatz |14| . and Tsimring and 
Pikovsky [H| studied the dynamics of a single noise ac- 
tivated bistable element with delayed feedback. Com- 
bining the properties of these two systems, Huber and 
Tsimring [2(| studied the properties of a globally coupled 
network of noisy bistable elements with uniform delays 
and derived a dichotomous mean field model based on 
the delay-differential master equation. Although for nu- 
merous systems the assumption of uniform time delays 
is justified [e.g. ElU^], most systems have time delays 
distributed over an interval rather than concentrated at 
a point [see l23i| . Thus, after a discussion of the dynam- 
ical properties of the network with uniform delays, the 
generalized case of distributed time delays is studied. 

This paper which is an extended version of |20| is orga- 
nized as follows: In the next Section the bistable-element- 
network is discussed for the case of uniform time delays. 
Two mean field models, namely the DZ-model (which for 
Gaussian processes reduces to the Gaussian approxima- 
tion) and the dichotomous model, are compared with the 
Langevin dynamics and their scopes of application are 
determined. The phenomenon of coherence resonance is 
discussed and a complete bifurcation analysis of the triv- 
ial equilibrium is carried out using a center manifold re- 
duction. In Sect. lIIII thc system dynamics is discussed for 
a discrete bimodal delay distribution. Then, in Sect. II Vl 
the model is further generalized, so that the mean field 
dynamics of a system with an arbitrary time delay distri- 
bution can be described. Finally, in Sect.^we introduce 
nonuniform coupling strengths which lead to new dynam- 
ical properties, such as amplitude death and chaos. 



II. UNIFORM DELAYS 
A. Langevin model 

The prototypical system considered here is modeled 
by a set of N Langevin equations, each describing the 
overdamped stochastically driven motion of a particle in 
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a bistable potential V = — x 2 /2 + x A /A, whose symmetry 
is distorted by the time delayed coupling to all network 
elements, 

JV 

x i (t) = x i {t)-x i {t) 3 + -^^2x j (t-T) + V2Dm, (1) 

3=1 

where r is the time delay, e is the coupling strength of 
the feedback and D denotes the variance of the Gaussian 
fluctuations which are mutually independent and 

uncorrected (&(*)&(*')) = tfytffjf - *')• 

The global coupling leads to an asymmetry of the local 
potential; that is, a positive feedback increases the prob- 
ability for an element to be in the potential well in which 
the majority of elements were at time t — r. The inverse 
holds for a negative feedback. 

System is explored numerically. In this paper, 
the numerical simulations are carried out using an Euler 
method. If not otherwise indicated the time-step and the 
number of elements are At = 0.01 — 0.05 and N — 2500. 

Our interest is mainly focused on the cooperative in- 
teractions of the individual network elements, i.e., on the 
dynamics of the mean field X = N^ 1 YljLi x j ■ F° r £ = 0, 
the elements are decoupled from each other. They jump 
from one potential well to the other randomly and in- 
dependently of each other. Therefore, in this case the 
mean field X = 0. For small |e|, the mean field remains 
zero. At a certain e = e st > which depends on the 
noise intensity D, but is independent of the time delay, 
the system undergoes a second order (continuous) phase 
transition and adopts a non-zero stationary mean field. 

For a negative feedback, a transition to a periodically 
oscillating mean field solution is observed at a certain e = 
£osc- < 0. Here and for the rest of this paper a (— /+) 
index means that the corresponding value is associated 
with a negative/positive feedback. 

Above a certain noise level Dh the transition at e sc- 
is second order as well. However, for D < Dh the sys- 
tem exhibits a first order (discontinuous) transition asso- 
ciated with hysteretic behavior. The noise intensity Dh 
depends on the time delay and is Dh = 0.07 for r = 100. 

For large time delays r 3> tk (jk is the inverse 
Kramers escape rate from one well into the other 0, 
l25jl. depending on the initial state the system adopts 
one of many accessible oscillatory states featuring dif- 
ferent periods. Even for a positive feedback, besides the 
stationary solution several oscillatory states with periods 
T < r are observed for e > e OS c+ ^ Sat- If the feedback 
is negative, the system only has oscillatory nontrivial so- 
lutions. The observed periods are T < 2r for e < £ OS c-- 

The basic dynamical states accessible by the system 
are illustrated in Figs.n an d|3 where the evolution of a 
single network element and the mean field are shown for 
different coupling strength. 



B. Gaussian approximation 

A mean field description for the dynamics of a glob- 
ally coupled set of thermally activated bistable elements 
with instantaneous interactions was proposed by Desai 
and Zwanzig |T3 |. For the sake of simplicity, we refer 
to this mean field description as the DZ-model. This 
model consists of a hierarchy of equations for the cu- 
mulants of the distribution function derived from the 
Fokker-Planck equation for the joint probability density 
function of all elements. Expressed in terms of moments 
M n {n = 1 . . . oo} the hierarchy assumes the simple form, 

M n =X(t- r)[4£»M„_ 2 + eM n _i] + M n - M n+2 , (2) 

where M_i = and Mq = 1. 

For large noise intensities, when the statistics of the 
individual elements are approximately Gaussian the hi- 
erarchy can be truncated (Gaussian approximation [see 
also|2(j). Applying this approach to our delayed feed- 
back system, the evolution of the mean field X is, in the 
Gaussian approximation, described by the following set 
of equations, 

X(t) = X(t)-X 3 (t)-3X(t)V(t)+eX(t-r), 

1 (3) 
-V{t) = V(t)-3X 2 (t)V(t)-W 2 {t) + D, 

where V = M 2 - Mf = N^ 1 J2?=i( x i ~ X Y is the vari- 
ance. 

To compare the theoretical predictions of the Gaus- 
sian approximation with the Langevin model Q we 
determine the maximum of the main peak in the power 
spectrum P pea k (see Fig. EJ. The evolution of the peak 
power as a function of the coupling strength can be used 
to study the Hopf bifurcation which describes the tran- 
sition to the oscillatory mean field regime. The pitch- 
fork bifurcation describing the transition to the station- 
ary mean field state is characterized by the dependence 
of the temporal mean of the mean field (X) t on the cou- 
pling strength. Fig. [3] shows the peak power P pea k and 
the temporal mean (X)t as a function of the coupling 
strength e for three different noise temperatures D. 

The phase diagrams of these models are shown in 
Fig. E] in the (D, e)-parameter plane. Fig. shows that 
away from the transition points the Gaussian approxima- 
tion correctly describes the Langevin dynamics. How- 
ever, near the bifurcation points the system dynamics 
is strongly non-Gaussian even for strong noise. Indeed, 
while the Gaussian approximation predicts that both bi- 
furcations are first order transitions (associated with hys- 
teretic behavior) over the entire noise range considered 
in this study (see Fig. 0J , the Langevin model produces 
first order transitions only for D < 0.7. 

The inclusion of higher order cumulant equations leads 
only to a slow convergence toward the true solution of the 
Langevin model. This is illustrated in Fig.[S] Thus, near 
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FIG. 1: Dynamics of a single network oscillator Xi and the network mean field X for different coupling strength e. The noise 
strength and the time delay is D = 0.1 and r = 100, respectively. For e < e OS c = —0.13 the mean field adopts a state of 
periodic oscillations. In the range e osc < e < e st — 0.11 the trivial equilibrium is stable. Finally, for e > e st the system adopts 
a non-zero stationary state. 
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FIG. 2: For e < 0, same as in Fig. Q but in the frequency domain. 



the transition points the DZ-model does not significantly 
simplify the Langevin description and the critical param- 
eters for the transition cannot be determined analytically. 



C. Dichotomous Model 

In order to describe the dynamics of the system near 
the bifurcation points we apply a dichotomous (i.e. two- 
state) approximation, which is complementary to the 
Gaussian approximation and which is valid in the limit 
of small noise, when the characteristic Kramers transi- 
tion time is tr- ^> 1, and small coupling strengths. The 
dichotomous theory neglects intra- well fluctuation of Xi. 



Thus, in the limit of small coupling, each bistable ele- 
ment can only take the values S1.2 = ±1- The collective 
dynamics of the entire network can then be described 
by the master equations for the occupation probabilities 
of these states m,2- This approach has been successfully 
used in studies of stochastic and coherence resonance [e.g. 
HI IT9IJ27I |2g| . For example, using this approach Jung 
et al. 0] i found nontrivial stationary mean field solutions 
in a globally coupled delay-free network of bistable ele- 
ments. 

The dynamics of a single element is determined by the 
hopping rates P12 and P21, i-e., by the probabilities to 
hop over the potential barrier from s\ to and from 
S2 to si, respectively. In a globally coupled system, ni.2 
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FIG. 3: The peak power P pea k and the temporal mean (X)t 
as a function of the coupling strength for the Langevin model 
(crosses), the Gaussian approximation (dashed line), where 
the double line indicates hysteretic behavior, and the dichoto- 
mous theory (solid line). The noise strengths is indicated 
in the upper right corner of each panel. The time delay is 
r = 100. For X = and D = 0.05,0.1,0.2 the Kramers 
times are tk = 659.4 , 54.1 , 15.5. 



and pi2,2i are identical for all elements and the master 
equations for the occupation probabilities read, 



"1 = -Pl2«-1 + P21™2 

n 2 = V\2n\ ~V2\ri2- 



(4) 
(5) 



The hopping probabilities pi2,2i are given by Kramers' 
transition rate |25j | for the instantaneous potential well, 



12,21 ^U"{ Xm )U"{x ) 

P K = ^ exp 



-AU 
D 



(6) 



where x m and xq are the positions of the potential min- 
ima and the top of the potential barrier, respectively. For 
our system, in the limit of small noise D and coupling 
strength e, they read [cf. 0], 



P12.21 



V2 =f 3 ai 
2vr 



■ exp 



1 =F 4«i 
AD 



(7) 



where a% = eX(t — r). 

As discussed above, the Langevin system either adopts 
a stationary or an oscillatory mean field state in the limit 
t — > co. Let us first consider the stationary case hip = 0. 
Making use of the probability conservation m + ri2 = 1, 
the occupational probabilities are given by 



ni.2 



P2142 
P12 +P21 



(8) 
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FIG. 4: Phase diagram for r = 100 of the Langevin model 
(crosses), the Gaussian approximation (dashed lines) and the 
dichotomous theory (solid lines and dotted lines). The solid 
line and the dotted line respectively depict the primary so- 
lution and the higher order solutions of Eq. I15H and 1161 . 
Phases separated by double lines indicate hysteretic behav- 
ior. For X — and D < 0.3 the Kramers time is tk > 10. 
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FIG. 5: The critical coupling of the Hopf bifurcation (e^sc < 
0) and the pitchfork bifurcation (e st > 0), respectively. Com- 
pared are the critical couplings resulting from the Langevin 
model (crosses) and the predictions of the DZ-theory (dia- 
monds) including different numbers of cumulants. For an 
even number of cumulants the DZ-theory predicts hystertic 
behavior which is not seen in the Langevin dynamics. 



Then, in the dichotomous approximation with si 2 = ±L 
the mean field X = sirti + S2^2 reads 



X = 7l2 — Til — 



Pl2 - P21 
Pl2 + P21 



(9) 



Substituting the hopping probabilities Q) into this 
equation yields the transcendental equation for the mean 
field magnitude 



X 



V2 - ZeX exp jeX/D) - V2 + 3sX exp (sX/D) 
V2 - SeX exp {eX/D) + V2 + 3eX exp {-eX/D) 

This equation always has a trivial solution X = 0, but 
for e > e a t it also has a pair of nontrivial solutions X = 



5 



±A It is easy to find A(e) for a fixed D numerically 
using IjlUI) . The critical value e st as a function of D for 
the pitchfork bifurcation, indicating the transition to a 
nontrivial stationary state, can be found analytically by 
expanding the r.h.s. of Eq. I|10(l at small X. This yields 
the following expression 



Est 



4L> 



3D 



(11) 



Let us now turn to the general case when X is al- 
lowed to be a function of time. Again, making use of the 
probability conservation and the expression for the di- 
chotomous mean field X = ri2 — ri\ we find the equation, 



X(t) =p\2 -P21 - (P21 +Pl2)X(t), 



(12) 



where the hopping probabilities P12.21 have the same 
functional form as in Eq. J2J, but now depend on the 
delayed mean field X(t — t) rather than X(t). 

To investigate the stability properties of the system, 
Eq. I|12|) is linearized about zero, 

X {t) = ^ exp(-l/4£>) ( e izJ£*(i - r) - X{t) 

(13) 

The characteristic equation for the complex eigenvalue A 
is found by making the ansatz X oc exp(Ai). It reads, 



£ (4-3J>) Ar _ 
AD 6 



(14) 



The trivial equilibrium loses its stability and undergoes 
a Hopf bifurcation indicating the transition to an oscilla- 
tory mean field state when the real part of the complex 
eigenvalue becomes positive. Therefore, the properties of 
the corresponding instabilities (i.e. frequencies and cou- 
pling strengths at the bifurcation points) can be found 
by substituting A = /i + iuj into Eq. 1)14(1 , separating real 
and imaginary parts and setting p = 0. This yields the 
following set of equations, 

v/2 

lot = cxp(— l/4:D)Tta,nuJT (15) 

(16) 



Est 



COS LOT 



This set of equations has a multiplicity of solutions, indi- 
cating that multi-stability occurs in the globally coupled 
system beyond a certain coupling strength. For finite 
time delays and positive coupling, besides the stationary 
solution, several oscillatory states with periods Tk close 
to but slightly larger than r/fc exist for e > £o SC+ {k = 
1,2,.. .}, where the transition points are ordered as fol- 



lows, < e st < e, 



OSC + 



< £osc+ ■ ■ • H the feedback is 
negative, the system has oscillatory solutions with pe- 
riods Ti close to but slightly larger than 2r/(2Z + 1) for 

P < £ l , , ., 

° ^ ^ — 1 1 " osc — ° osc - 

00, the transi- 

-0 _ 4D 



{I = 0,1,...}, where > e° sc _ > e, 
In the limit of large time delays r 



with the corresponding periods being Tk — > r/fc and 
Ti -> 2t/(2Z + 1), respectively. 

In order to compare the predictions of the dichotomous 
model with the Langevin dynamics, the peak power P pea k 
and the temporal mean (X) t , resulting from the dichoto- 
mous theory, are also plotted in Fig. The phase dia- 
gram for the dichotomous theory is shown in Fig. 0] 

Fig. |3 and 01 show that the dichotomous theory agrees 
with the Langevin dynamics quite well for small noise in 
the range D « 0.07—0.3 in the neighborhood of the bifur- 
cation points. The theory also correctly describes the bi- 
furcation type. Indeed, the dichotomous theory predicts 
accurately the noise strength £>h (= 0.07 for r = 100) 
at which the Hopf bifurcation changes from supercriti- 
cal (second order) to subcritical (first order). However, 
for very small D the Kramers time becomes very large, 
and the accuracy of numerics becomes insufficient for a 
comparison with the theory. 



D. Complete bifurcation analysis 

A complete bifurcation analysis of the trivial solution 
X = of Eq. i|13[l in the (D, e, r)-parameter space can be 
accomplishe d by carrying out a center manifold reduction 
[see e.g. |2|| l30| : that is, the normal form coefficients of 
the bifurcations in our dichotomous mean field model can 
be expressed in terms of the system parameters. 

For a general class of delay differential equations of the 
form, 



x(t) 



x(t) + -fix(t - 
+^lix{t)x{t - 



7 3 cc(t) x(t 
^3 



■r) 
(17) 



such a reduction to normal forms of the pitchfork and 
Hopf bifurcations has been carried out in Refs. [sil l32l | . 
If we cast the equation for the mean field dynamics of 
our model in this form, we can use the results in Refs. 
[3TL to determine the functional dependence of the 
normal form coefficients on the parameters D, e, and r. 
This can be achieved by a series expansion of Eq. i|12|) 
up to the third order and rescaling of time. 

The normal form of the pitchfork bifurcation reads 



i = az + bz , 

where z is a coordinate on the center manifold, 
normal form coefficients are 

£ ~ Est 



b = 



£ st (l 
Bi- 



-To) 

12DB 2 



384(1 -t )L> 3 ' 



(18) 
The 

(19) 
(20) 



tion points s. 



OSC+ 



e st and 4 



4-3Z5 



where B x = e 3 (81L> 3 + 108Z? 2 + 144D - 64), B 2 = 
£ 2 (9-D 2 + 24L>-16) and r = -V2cxp(-l/4L»)r/7r. Set- 
ting a = and solving Eq. (|19|l for e we again find the 
critical coupling of Eq. I|ll|) . One can show that b < 
for D > and e = e s t (i.e. a — 0). Consequently, the 
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FIG. 6: Stability diagram for the trivial equilibrium resulting 
from the analysis of the pitchfork bifurcation. The solid and 
the dashed lines respectively depict the b=0 and a=0 contour 
lines. The stability diagram for the pitchfork bifurcation is 
time delay independent. 



pitchfork bifurcation at e st is always supercritical. The 
stability diagram resulting from center manifold reduc- 
tion for the pitchfork bifurcation is shown in Fig. 

The normal form of the Hopf bifurcation in polar co- 
ordinates r and 9 on the center manifold reads 



ar 



pr 



(21) 



The coefficients determining the stability of the trivial 
equilibrium and the order of the Hopf bifurcation are /i 
and a [3j|. Expressed in system parameters, they read 



/* = — 



1 



1-TQ 



sm tp 



a 



£ st cos ip J \(1 
5x53-52(1-370 



- 2 cos 2 ip) 



128([l-r 



(COS 2 if - 



<p 2 )D 2 
t )/(4D cos ip 



,(22) 
(23) 



tq and ujq 



where B3 
tan ip. 

Setting fi = and solving Eq. (1221 for e yields the crit- 
ical coupling as a function of the noise strength e osc (D) 
which coincides with Eq. I|16|) ■ Setting the first Lyapunov 
coefficient a = 0, we can find £ a =a(D). The two func- 
tions £ QSC {D) and e a= o(D) intersect at a noise level 5h 
denoting the parameter values for which the Hopf bi- 
furcation changes from supercritical to subcritical. The 
stability diagram resulting from the analysis of the Hopf 
bifurcation is shown in Fig. \7\ 

Let us now discuss the bifurcation properties in the 
limit of large and small time delays as well as vanishing 
noise and compare them with those of a single oscillator 
system. The critical coupling e st of the pitchfork bifurca- 
tion is time delay independent and goes to zero for van- 
ishing noise. However, the critical coupling of the Hopf 
bifurcation depends on the time delay (see lower panel in 
Fig. UJ. As the time delay increases, the maximum of the 
primary Hopf bifurcation line £o SC _ approaches the origin 
in the (e, D) plane meaning that oscillations may occur 
at an arbitrary small feedback strength for the properly 
tuned noise level. This should be contrasted to the dy- 
namics of a single noise-free oscillator with time-delayed 
feedback that only exhibits oscillations at strong negative 
feedback (e < —1). For very small time delays t — > 0, 
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FIG. 7: Primary solutions of fx = (Eq. 133) and a = 
(Eq. I23H in the (s, D)-parameter-space. Upper panel: The 
boundaries /i = (dashed line) and a — (solid line) for 
a system with r = 100. The black dashed line and the gray 
dashed line depicts the parameter values of the subcritical and 
supercritical Hopf bifurcation, respectively. The lower panel 
shows the same curves for a system with r = 10 (dashed line), 
r = 100 (black solid line) and r = 1000 (gray solid line). 



E. Coherence resonance and system size effects 

The system studied in this paper exhibits the phe- 
nomenon of coherence resonance [e.g. I3"H l3fl l3fl l37| and 
array-enhanced resonance [2fl l3Sj . 

Let us discuss this in turn. If our system adopts an os- 
cillatory state, the double-well potentials of the elements 
are tilted asymmetrically, due to their coupling to the de- 
layed mean field; that is, the potential barriers separating 
the two wells are periodically rising the lowering. If the 
period of this oscillation T matches the time scale tk of 
the noise- induced inter- well fluctuation, i.e., if the mean 
field oscillations synchronize with the hopping rate, we 
can expect that the number of elements contributing to 
the oscillation and consequently the order of the oscilla- 
tory state reach a maximum. In this spirit the time scale 
matching condition for such a synchronization, which is 
given through 



2t k = T, 



(24) 



=Foo. 



is a reasonable condition for the maximum order of the 
oscillatory state |28j . 

To quantify the order (i.e. coherence) of the oscilla- 
tory state we introduce the coherence parameter (3 = 
Hujpeak/ Alu, where H is the height of the main spec- 
tral peak at w p0 ak and Aw is its halfwidth. Using the 
Langevin model l[T)l. the coherence measure (3 is deter- 
mined as a function of the noise strength and in Fig. [S] 
compared for systems of different size N . 

Clearly, the coherence curves have a maximum. The 
noise strength maximizing the coherence is D$ w 0.08. 
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FIG. 8: The coherence of the oscillatory states (3 as a func- 
tion of the noise strength D for systems of different size N. 
The time delay and the coupling strength are r = 100 and 
e = —0.2, respectively. Left panel: The coherence of the mean 
field oscillations. Right panel: The coherence of a single ele- 
ment Xi out of the N network elements. 
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FIG. 9: The Kramers time tk and the half the period of 
the mean field oscillations (in units of the time delay) as a 
function of the noise strength. The parameters are r = 100 
and e = -0.2. 



This noise strength can also be derived from the time 
scale matching condition in Eq. (|24|l . The Kramers time 
tk = 1/p is given through Eq. Q and the period of the 
oscillations T beyond the critical coupling can be deter- 
mined numerically. In Fig.[5]the two time scales are plot- 
ted as function of the noise strength. The curves intersect 
at D=0.08 substantiating the consistency of theory and 
Langevin model. 

The resonance curves in Fig. [5] show that the coher- 
ence of the oscillatory states increases with increasing N, 
a property which was reported for other systems and is 
sometimes referred to as array- enhanced resonance |38j . 
Interestingly, the enhancement of the temporal regularity 
with increasing system size is only observed for macro- 
scopic mean field oscillations, while the inverse holds for 
"subcritical coherence" . That is, the coherence observed 
in the power spectra of subcritical mean field fluctuations 
(i.e., for |e| < |e OS c±|) decays inversely proportional to the 
number of elements in the network, and becomes negligi- 
ble for N > 10. This is show in Fig.^JI Qualitatively, the 
same dependency on the system size is found if the de- 
layed average does not include the delayed element itself, 
i.e., the element xi is coupled to Xi(t— r) = YljLi^j^ti x j- 
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FIG. 10: The coherence /3 as function of the coupling 
strength. For e < and e > the spectral peak frequency is 
/peak = ^ pC ak/27T w 0.5 1/t and / pcak ?s 1.0 1/r, respectively. 
The dash-dotted vertical line depicts e^sc- an d consequently 
separates domain of the macroscopic (i.e. supercritical) mean 
field oscillations from the domain of subcritical coherence. 
The right panel shows the same as the left, but has a logarith- 
mic scale for /3, which helps to uncover the weak subcritical 
coherence properties. 



III. TWO DELAYS 

A. Langevin model 

We want to generalize the above system by introduc- 
ing multiple time delays and non-uniform coupling terms. 
Let us carry out the generalization progressively and 
study first the dynamics of a bistable element network 
with a discrete bimodal delay distribution (i.e. with two 
time delays) and uniform coupling. Assuming that the 
time delay of the interaction between two elements is 
entirely determined by the "transmitting" element the 
system dynamics is described by the following set of 
Langevin equations, 

Xi - Xi-x?+^Xx(t-Tx) + ^X 2 (t-T 2 ) + V2D£(t) (25) 



where 



and 



9 N/2 



2 N 
X ^=N £ x i®> 

j=N/2+l 



(26) 



(27) 



are the mean fields of the elements associated with time 
delay n and T2, respectively. Here, it is assumed that 
the number of oscillators is the same in both group. 



B. Dichotomous theory 

We want to use the dichotomous theory in order to 
study the mean field dynamics of model (|25|l . Thus, the 
theory developed in Sect. Ill Cl has to be extended accord- 
ingly. In order to describe the collective dynamics of the 
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two-delay system, two equations are needed, respectively 
describing the mean field evolution of the oscillator group 
associated with n and r 2 , 



Xi t2 (t) = P12 ~ P21 ~ (P12 + P2l)^l,2(*)< 



(28) 



The mean field of the entire system then is, X — (X\ 
X 2 )/2, and the hopping probabilities are given by, 



P12.21 



V2 jjj 3 « 2 
2tt 



■ exp 



1 T 4a 2 
4D 



(29) 



where a 2 = e [X x {t - t\) + X 2 (t - r 2 )] /2. As for the 
model with a single (i.e. uniform) time delay, the nu- 
merical integration of the Langevin system l|25l) reveals 
pitchfork and Hopf bifurcations describing the transitions 
to nontrivial stationary states and oscillatory states, re- 
spectively. The critical couplings for the bifurcations can 
be found with a linear stability analysis of Eqs. (|28|l near 
the trivial state X = 0. The procedure, which is analo- 
gous to the stability analysis carried out in Sect. Ijll Cfl . 
yields the transcendent equation for the complex eigen- 
value A, 



A = 



S ± VS 2 - 4A 



(30) 



Here, S and A respectively are the trace and the deter- 
minant of the Jacobian matrix 



J = c 



gi + d g 2 

gi g 2 + d 



(31) 



where the matrix elements are given through c = 
-V2exp(-l/4L>)/87rL>, g h2 = e(3D - 4) exp(-Ar lj2 ) 
and d = 8D. For a positive coupling Eq. (|3(JI) has al- 
ways a real eigenvalue. At a certain critical coupling 



AD 



Est 



4-3L> 



(32) 



the eigenvalue becomes positive indicating the pitchfork 
bifurcation. This bifurcation is time delay independent 
and is thus identical with those found for the system with 
uniform time delays [cf. Eq. 

For finite f = (n + r 2 )/2 and e, Eq. (|3Tj|l possesses 
also a finite number of complex solutions. The critical 
couplings of the corresponding unstable modes (i.e. of 
the Hopf bifurcation) are given by the following set of 
equations, 



V2 



lot = exp(— 1/AD)t tanwr, 

7T 



8DTTLU 



(3D - 4)(v / 2exp[-l/4L>] J s - itljJ c ) 



(33) 
(34) 



Here 



J s = — (sin ujti + sinwr 2 ) = sin lot cos lo a (35) 
J c = —(cos ojti + cos ujt 2 ) = cos lot cos uia, (36) 
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FIG. 11: Phase diagram of the globally coupled two-delay 
network determined using the dichotomous model (solid lines) 
and numerical simulations of the Langevin model (markers). 
The phase diagram is show for different a = |ri — T2|/2. The 
mean time delay is f = 100. 



where a = \t\ — r 2 \/2. The above set of equations for the 
critical coupling is the two-delay analog to Eq. (|15fl and 
(|16fl . Again, we find a multiplicity of solutions, leading 
to the multi-stability of the system in a certain area of 
the parameter space. Furthermore, Eqs. H33|) . l|34| show 
that while the frequencies of the oscillatory states only 
depend on the mean time delay f , the critical coupling 
strengths of the Hopf bifurcations depend additionally on 
a. 



C. Phase diagrams 



The phase diagram and the frequencies of the unsta- 
ble oscillatory modes of the two-delay system are theo- 
retically determined using Eqs. I|32|l - I|34|l and compared 
with numerical findings resulting from simulations of the 
Langevin model l|25|l . The phase diagram is shown in 
Fig. 1111 for different a . The phase diagrams including 
higher order solutions of Eq. (|33|) and (|34|l are presented 
in Fig. E| Also, the frequencies of the corresponding 
unstable modes (which are er-independent) are shown 
in this Figure. The Figures show that near the bifur- 
cation points the predictions by the dichotomous the- 
ory are reasonably good for weak noise in the range 
(0.07 < D < 0.3). 

Furthermore, we find that the first bifurcation of the 
trivial equilibrium at e > is always a pitchfork bifurca- 
tion. The first bifurcation at e < is a Hopf bifurcation, 
which for er < 30 is determined by the primary solution 
of Eq. l(33|) and while for a > 30, depending on 

the noise intensity, the first transition may also be deter- 
mined by higher order solutions associated with higher 
frequencies. 
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FIG. 12: Upper panel and two lower left panels: Phase diagrams of our globally coupled two-delay network with f = 100 and 
a = 0, 10, 20, 30, 40. The green line depicts the critical coupling of the pitchfork bifurcation and the other lines depict those 
of the primary Hopf bifurcation as well as some higher order solutions (i.e. solutions 1-15) of Eq. 1331 and l3~4l The markers 
depict the first bifurcation at e < and e > resulting from numerical simulations of the Langevin model (in these simulations, 
starting with e = 0, the coupling strength is increased until a bifurcation occurs). Matching colors of markers and lines mean 
that the bifurcation type and associated frequency are in agreement. Lower right panel: The frequencies of the corresponding 
unstable modes. They do not depend on a, but slightly vary with the noise strength D. 



IV. MULTIPLE DELAYS 
A. Langevin model 

In this Section we further generalize our delayed feed- 
back system by introducing multiple time delays and 
study the stability properties in dependence of the sta- 
tistical moments of an arbitrary time delay distribution. 

The general Langevin model with many time delays 
reads 

N 

x t = Xl - + x 3 (* ~ T v ) + (*) • ( 3? ) 

3=1 

Such general models in which the time delays depend 
on both the "transmitting" and the "receiving" element 
cannot directly be described in terms of a mean field 
theory. However, the system becomes mathematically 
tractable if we assume that the time delays only depend 
on the transmitting elements j, 

N 

±i = x t - xf + j- x 3^ ~ t j) + V^5$(t). (38) 

3=1 

In order to check if such a simplification is justified, 
numerical simulations of model 13 7|) and i)38|) are carried 
out and compared. In these simulations the distribution 
of the time delays is Gaussian, i.e., it is fully determined 
by its mean f and standard deviation a. Fig. 1131 compar- 
ing the critical coupling strength of the Hopf bifurcation 
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FIG. 13: The critical coupling of the Hopf bifurcation as a 
function of the noise strength D for different a of the Gaussian 
time delay distribution with r = 100. The markers and the 
solid lines depict the critical couplings resulting from model 
13711 and model 1381 . respectively. 



for different a, suggests that the above simplification is 
justified in order to study the stability properties of a 
bistable-element-network with time delays. 

This surprising result not only renders possible an an- 
alytical description of networks with distributed delays 
but also implies that the number of operations, which 
have to be carried out to study such systems numerically, 
can be reduced from 0{N 2 ) to O(N). 
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B. Dichotomous theory 

Let us now develop the dichotomous theory for 
the globally coupled bistable-element-network with dis- 
tributed delays. 

For that purpose we coarse grain system i|38[l . The 
coarse graining is accomplished as follows: The range of 
possible time delays is divided up in M intervals I k {k — 
1,2, ... , M}. The size of the intervals A k is chosen, so 
that the number of bistable oscillators associated with a 
delay, fitting in a particular interval, is for each interval 
the same m — N/M. In this way oscillator groups are 
formed whose mean field can be expressed as, 



m — ' 



(39) 



where I k = [r k , r k+1 [, r k = J2i=i A * and J = 1 ■ ■ • N - 

Assuming that -C f/a, where f and a are the mean 
and the standard deviation of the time delay distribution, 
Eq. i|38|) can then be approximated by, 

M 

x i = x i -xf + -^J2 n k{t-T k ) + V2Dat)- (40) 
fc=i 

The master equations expressing the dynamics of sys- 
tem 1)400 in terms of occupation probabilities read, 



n 2 M = Pi2»i,fc ~P2iri2,k- 
Here the hopping probabilities are given by, 
n/2 T 3a 3 / 1 T 4a 3 



Pl2,21 



2tt 



where a 3 = (e/M) Y,k=i ^k{t - r k ). 

For large oscillator groups (m — > oo), fi^ = ni^si + 
^2,fcS2 = «2.fe — ?ii,fc holds. With this and the probability 
conservation ni^ + ri2,fc = 1 we can find the following set 
of equations: 

fik(i) =pi2 -P21 - (P21 +Pi2)O fc (i). (44) 
The Jacobian matrix of this system is given through, 



■ exp 



AD 



(41) 
(42) 



(43) 



f9l +d g 2 
9i 92 + d 



9m 
9m 



(45) 



9i 



92 



9M 



where c = -V2exp(-l/4D)/(4M7rD), g k = e(3D - 
4)exp(— Xr k ) and d = AMD. With this Jacobian the 
characteristic equation, determining the stability of the 
trivial equilibrium X — 0, becomes 



(dc - A) 



A/-1 



k=l 



(46) 



Setting A = and solving Eq. I|46[) for e yields the 
critical coupling for the pitchfork instability, 



Est 



AD 
4- 3D' 



(47) 



It is time delay independent and thus identical with that 
found in previous Sections of this paper. 

The properties of the Hopf bifurcation (i.e. the fre- 
quencies of the unstable modes u and the critical cou- 
pling e osc can be found by substituting A = /i + \ui into 
Eq. (|46[) . separating real and imaginary parts and setting 
/i = 0. This yields, 



— exp(-l/4D)f^ 

7T L 



where 



h = 



^ M M 

— > sin ujT k , Ic — — > 



COS UJT k . 



(48) 



(49) 



k=l 



k=l 



For large systems N — > oo, the number of groups M — > oo 
and thus 

oo oo 

Z s = y P(t) sin Lordr, I c = J P(t) cos curdr, (50) 



where P(r) is the time delay distribution function. 

We can express the time delay distribution function in 
terms of cumulants K n (3^, |4jj and solve the integrals in 



I s = sin(gi) exp(# 2 ), / c = cos(gi) exp(g 2 ), (51) 



where 



.91 



.92 



2m+l 



Consequently, 



'-^ i(2m + 1)! 

m=0 V ; 
m=l v ; 



tan(gi). 



2m+lj 



(52) 
(53) 

(54) 



Since for symmetric distribution functions all odd cu- 
mulant moments except the first one K\ = f are zero, 
I s /I c = tanwr holds. That is, in the case of a symmetric 
distribution of the time delays, the frequencies of the un- 
stable modes in Eq. I|48l) depend only on the mean time 
delay. 

Let us now determine the critical coupling of the Hopf 
bifurcation. For large time delays r»rj( the low-order 
solutions of the transcendental equation l|48|l yield fre- 
quencies lu <C 1. Thus the real part of equation 146|) can 
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FIG. 14: Phase diagram of the globally coupled bistable- 
element-network with uniformly distributed time delays de- 
rived from the theoretical model (solid lines) and numerical 
simulations of the Langevin model (markers). The phase dia- 
gram is shown for different standard deviations a of the delay 
distribution function. The mean time delay is f = 100. 



be linearized near uj = and the critical coupling of the 
Hopf bifurcation becomes, 

£ ° SC = (3D - 4) (£V2exp(-l/4£>)J 8 - [l - J] W c ) ' 

(55) 

Then, for large systems N — * oo the critical coupling is, 



4L> 



(4 - 3D)I C ' 



(56) 



with I c = 3 sm(u;T) sin(5wer/3)/(5aj<7) and I c — 
cos(wf ) exp(— w 2 cr 2 /2) for uniform and Gaussian distri- 
butions, respectively. 



C. Phase diagrams 

Eqs. and iJSBJ are used to determine the 

phase diagram and the frequencies of the unstable oscil- 
latory modes / = uj/(2tt) of a bistable element network 
with uniformly distributed time delays |47| . The theoret- 
ical predictions are compared with numerical simulations 
of the Langevin model l|37l) . The number of bistable el- 
ements in these simulations is TV = 300. The results are 
shown in Figs. 1141 and!15l 

Again, we find that near the transition points and for 
weak noise intensities the predictions of the dichotomous 
theory are reasonably good. Consequently, the Langevin 
models (|37|) and l|38|l are in this regime equivalent as 
regards the dynamical properties of the mean field. 

Eqs. I|48l) and (|56[) have a multiplicity of solutions 
meaning that multistability is also present in our sys- 
tem in the limit of continuously distributed delays. The 
bifurcation diagrams including the higher order solutions 
are shown in Fig. 1161 The Figure shows that unlike the 
two-delay system, the first transition at e < is always 
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FIG. 15: The frequencies of the unstable modes at the bifur- 
cation points resulting from the Langevin model (markers) 
and the the dichotomous model (solid line), which are (see 
Eq. 1481 independent of a. For uniform and Gaussian distri- 
butions the frequencies depend only on the mean time delay 
(see Ea.l48land l54l . 
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FIG. 16: Same as in Fig. 1121 but this time for networks 
with uniformly distributed time delays with r = 100 and 
cr = 0, 10, 40 



determined by the primary solution associated with the 
frequency / w 0.5 r. 

The comparison of the phase diagrams for delay dis- 
tribution functions of different widths a shows that the 
regions of oscillatory states in the parameter space are 
reduced with increasing a . This trend was already appar- 
ent in the two-delay system, although less pronounced. 
These findings suggest that nonuniformity of the time 
delays inhibits the occurrence of Hopf bifurcations and 
consequently increases the stability of the trivial equilib- 
rium. 

Eventually, we like to mention that the coherence reso- 
nance phenomenon discussed in Sect. UTEl is also present 
in systems with multiple delays in the oscillatory domain 
of the phase diagram. 
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NONUNIFORM COUPLING 



A. Langevin model 



The collective dynamics of the bistable-element- 
networks described above, is restricted to periodic os- 
cillations and stationary states. In this Section we want 
to check whether the complexity of the dynamics is in- 
creased if instead of the uniform coupling, non-uniform 
couplings are applied. To this end, we extend the two 
time delay model (|25|l by introducing two different cou- 
pling strengths. The Langevin equations of the new 
model read, 



e -l Xl (t-n) + ^X 2 (t-r 2 ) 



(2 



Xi(t-n) 



X 2 {t 




where i = 1, . 
mcnts Xj and 



2D£(t), 
(57) 

, N/2 and j = N/2 + l,...,N. The ele- 
belong to a group of bistable oscillators 
which are associated with time delay t± and t 2i respec- 
tively. It is assumed that the two groups are of equal 
size. The above set of equations describes a system in 
which each element couples to all the elements belong- 
ing to the same group with a coupling strength E\ and 
to all the elements of the other group with e 2 ; that is, 
the two coupling parameters indicate the strength of the 
intra-group coupling (ei) and inter-group coupling (£2), 
respectively. 

B. Dichotomous model 

We apply the dichotomous theory to system l(57|) and 
proceed in a manner analogous to the previous sections. 

The evolution of the mean field of each group of oscil- 
lators is described by 

Xi, 2 {t) =p\ 2 2 -P21 2 - (P12 +P 1 2i ) x i-2( t )- (58) 
Here, the hopping probabilities are, 



Pl2,21 



-Pl2,21 



V2 T 3a 4 

2tt 
V2 T 3a 5 



■ exp 



■ exp 



1 =F 4«4 \ 
AD J 



(59) 
(60) 



2tt "V 4L> 
where 

a 4 ,5 = [siX h2 (t - n ) + s 2 X 2tl (t - t 2 )] /2. (61) 

Next a linear stability analysis is carried out. The lin- 
earization of Eq. 1581 about the trivial equilibrium yields 
the Jacobian, 



J 



([3D - 4]ei 



d [3D-4]e 2 e" Ar2 , 



I [3D - 4]e 2 c" Ari [3D - 4]eie" Ar2 



where c and d are the same as in Eq. 1 13 It . 

Substituting the trace S and determinant A of the Ja- 
cobian matrix 1)62(1 into the characteristic equation 



A 



S ± y/S 2 - 4A 



(63) 



and keeping the intra-group coupling strength E\ fixed, 
yields the critical coupling for the pitchfork bifurcation, 
which can occur for positive and negative inter-group 
feedbacks, 



st 
e 2 



± 61 



8D 



3D -4 



(64) 



In order to find the critical values for the Hopf bifurcation 
E2 SC , we substitute A = fj, + iu> into the characteristic 
equation and set /i = 0. Then, the separation of real and 
imaginary part yields the two equations, f r (uj,e 2 ) = 
and fi(u>, e 2 ) — 0, where 

E 2 

f r (u,e 2 ) = E 1 e 1 J s Lu + ^-(e 2 1 -el)cos(2ujf) 



+ E 2 £ 1 J c + ^eM-l/2D)-4u 2 , (65) 
E 2 

fi{cj,e 2 ) = E l e 1 J c w+ -j- (4 - £2) sin(2wf ) 

+ ^ 2 £iJ s + ^^exp(-l/4Z?). (66) 



Here, E 1 = V2(3D - 4) exp(-l/4L>)/(vrL>) and E 2 = 
i?iv / 2exp(— 1/4D)/tt. The terms J s and J c are given by 
Eq. pB*|l and respectively. 

For finite f and £ 2 the above set of equations has a 
finite number of roots (£ 2 sc ,w), which can be found nu- 



merically. 



C. Phase diagram 



(62) 



In order to explore the dynamics of the system with 
two coupling strengths we carry out numerical simula- 
tions of the Langevin model l|57(l and compare the results 
with the theoretical predictions derived in the previous 
Section. In these simulations the strength of the intra- 
group coupling £1 and noise D are chosen so that in the 
absence of inter-group couplings £ 2 = the mean fields 
of the two oscillator groups X\ and X 2 oscillate inde- 
pendently with frequencies /1 ~ 1/2ti and f 2 w 1/2t2 
respectively (cf. Sect. Ill CI . We may then expect that for 
I £2 1 > the system reveals dynamical properties reminis- 
cent of those of two coupled, limit-cycle oscillators, such 
as chaotic behavior 1 121 143| and the amplitude death 
phenomenon || |2^, |44| . Fig. El shows time delay rep- 
resentations of the time series of Xi(t) for inter-group 
couplings of different strength £2. In certain regions the 
system indeed shows the amplitude death phenomenon 
and in the range < |e| < £ c haos irregular motions are 
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observed. Numerical evidence suggests that these mo- 
tions are chaotic. Indeed, the time series analysis yields 
broad band power spectra and positive maximum Lya- 
punov exponents. The determination of the Lyapunov 
exponents is below discussed in greater detail. 

The comparison of the phase portraits in Fig. 1 171 shows 
slight deviations between theoretical predictions and the 
Langevin dynamics (e.g. for £2 = —0.1). These devi- 
ations stem from the elimination of the noise fluctua- 
tions in the dichotomous model and different phase shifts 
between the two oscillator groups X\ and A 2 . How- 
ever, the predictive power of our model is confirmed in 
Fig. [^J where the theoretical peak power P pea k and the 
corresponding period T pea k in dependence of the cou- 
pling strength ej are compared with those resulting from 
Langevin simulations. 

Let us now explore the phase space of the system with 
nonuniform couplings in greater detail. 

The phase space regions of nontrivial stationary states 
are determined by Eq. and those where mean field 
oscillations and amplitude death occur, are given by the 
roots of Eq. I|t)5|) and (|66[) . These roots are determined 
numerically. We find that the solutions of Eq. 1|33[1 are a 
subset of the the solutions of / r ,i(w) = 0. Thus, the cor- 
responding critical values mark boundaries which quali- 
tatively arc equivalent to those found in previous models. 

However, Eq. I165H also yields new solutions marking 
the boundaries between the the zones of amplitude death 
and the areas of nontrivial dynamics in the presence of 
weak inter-group couplings. Within this areas there may 
occur islands of chaotic dynamics. Indeed, an analysis of 
the mean field evolution yields positive Lyapunov expo- 
nents for < |e| < £ c hao S - 

Since intrinsically our time delay system is infinite di- 
mensional the maximum Lyapunov exponents are here 
determined by an analysis of the time series resulting 
from Eq. I|58|) . The analysis is carried out using tools 
provided by the TISEAN software package [ill , liq . 

As stated above this process yields in some phase space 
regions clear evidence of positive maximum Lyapunov 
exponents in the range < A [1/f] < 0.03. 

The phase diagrams illustrating the different dynamic 
regions are shown in Fig. ^| The Figure shows that 
chaotic dynamics only occurs for strong intra-group cou- 
plings £1 > 0.4, i.e., when the individual oscillations of 
the two groups are strong enough. 



a variety of stable oscillatory mean field states are ac- 
cessible for a positive as well as negative feedback. The 
coherence of the oscillatory states is maximal for a certain 
noise level; i.e., the systems demonstrate the coherence 
resonance phenomenon. 

For symmetric time delay distributions the frequencies 
of the oscillations depend only on the mean time delay. 
However, the critical couplings of the corresponding Hopf 
bifurcations depend also on the higher order cumulants of 
the time delay distributions. Indeed, our findings suggest 
that nonuniformity of the time delays inhibits the occur- 
rence of the Hopf bifurcations and consequently increases 
the stability of the trivial equilibrium. This may be im- 
portant for time delay systems such as neural networks 
and genetic regulatory networks, since the degree of time 
delay nonuniformity, which is often related to the diver- 
sity in the connectivity of the underlying network, affects 
the accessibility of the nontrivial dynamical states. 

The dichotomous theory based on delay-differential 
master-equations, which has been developed in this ar- 
ticle, adequately describes the bifurcations of the triv- 
ial equilibrium in the limit of small noise and coupling 
strength. Furthermore, the theory allows for the appli- 
cation of a center manifold reduction and thus for a com- 
plete bifurcation analysis of the trivial equilibrium. Far 
away from the bifurcation points the mean field prop- 
erties are well described by a Gaussian approximation. 
However, a theoretical approach for the description of 
the dynamics in the regime of strong noise near the tran- 
sition points is still lacking. 

The collective dynamics of the networks of bistable el- 
ements with uniform coupling strength is restricted to 
periodic oscillations and stationary states. However, our 
model with nonuniform coupling strengths shows that for 
certain coupling distributions, the system behaves like 
a network of coupled limit cycle oscillators and conse- 
quently, demonstrates in certain parameter-space areas 
the amplitude death phenomenon or exhibits a chaotic 
evolution of the mean field. 

This paper discusses the dynamics of globally coupled 
systems with time delays. However, in many systems the 
connectivity is sparse. Since this is a particular case of 
systems with nonuniform coupling we may expect that 
this endows the system with more complex dynamical 
properties. This issue should be addressed in future stud- 
ies. 



VI. SUMMARY AND CONCLUSIONS 

The dynamics of networks of noisy bistable elements 
with time-delayed couplings was studied analytically and 
numerically. Depending on the noise level, the systems 
undergo ordering transitions and demonstrate multista- 
bility; that is, for a strong enough positive feedback the 
systems adopt a nonzero stationary mean field state, and 
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FIG. 17: Time delay representations (phase portraits) of the evolution of the mean field Xi in dependence of the inter-group 
coupling strength £2- Shown are the evolutions resulting from the Langevin model (upper two rows, Eq. I57H and the mean 
field model (lower two rows, Eq. I58H . The parameters are: D — 0.1, n = 60, Ti = 140, Ei = —0.4. 




FIG. 18: The peak power P pca k (upper row) and the cor- 
responding period T pca k = 2n/u pcs _] i (lower row) of the two 
oscillator groups X\ and X2 resulting from simulations of the 
theoretical mean field model (I58H and the Langevin model 
1571 1 . respectively. The parameters are the same as in Fig. 1171 
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